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Abstract. The existing literature on stochastic simulation of chemical reaction networks has a tendency to move as quickly as possible to the 
abstract formulation of the stochastic dynamics in terms of probabilities based on the concept of the Chemical Master Equation (CME), largely 
ignoring sample path representation. In this publication we discuss both theoretical basis and numerical approach for the problems in this area 
using sample path methods as a crucial part of the process. Relying as it does on a representation of the underlying stochastic processes as a weak 
solution of a system of stochastic differential equations driven by Poisson random measures this approach brings to bear a heretofore ignored but 
quite effective problem solving methodology. 

We first present a simple and intuitive way of partitioning species and reactions of the interaction network into different groups. We then discuss 
how original stochastic dynamics with state dependent intensities of transitions can be reformulated in terms of jump-diffusion stochastic differential 
equations driven by both Wiener noise sources and Poisson random measures. Finally, we show that this approach facilitates the construction of 
hybrid simulation techniques, an important step in the creation of efficient techniques for modeling multi-scale stochastic dynamics of the reaction 
networks. Numerical methods related to sampling events from Poisson random measures are demonstrated on simple intuitive examples. Error 
control analysis of the finite differences scheme is also presented. 

Key words. Stochastic algorithms, chemical reaction networks, Ito-Skorohod stochastic differential equations, Poisson random measure, 
hybrid systems 
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1. Introduction and Motivation. Many important problems in study of complex systems require consideration 
of dense and nonlinear interaction between different functional modules. Hybrid system is a framework first proposed 
by the engineering community |Sas99| to model systems that exhibit both continuous and discrete state changes, and 
has, in recent years, found increasingly wide applications in many practical fields in connection to both complex man- 
made systems or physical systems found in nature. However, most of the models proposed in literature so far are 
deterministic, and are not suitable to model inherent randomness. 

New data from high-throughput biology is enabling more ambitious, complete and validated simulation biological 
models and a lot of attention has recently been paid to the development of the computational methods for the in 
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silico investigation of complex biological processes on a cellular level. To be valuable to biological and biomedical 
research computational methods must be scaled to account for to large number of interacting component of different 
nature and complex pathways. As it was already noted by many researches stochastic effects play an important 
role in the functionality of such systems | ARM98 1 and approximation of molecular behavior by continuous models 
fails to reproduce many interesting biological phenomenal RWA02|. The deterministic approach can be particularly 
misleading if the molecular population of some critical reactant species becomes so small that microscopic fluctuations 
can conspire with reaction channel feedback loops to produce macroscopic effects. Recent work shows that this can 
happen with dramatic consequences in the genetic/enzymatic reactions that go on inside a living cell | RWA02 1 or 
non-equilibrium combustion reactions. 

Correct methods for performing stochastic simulation of the discrete event simulation are usually based on 
the Kinetic Monte Carlo Methods(KMC) IBML75l.Stochastic Simulation AlgorithmfSSAI IGii??! lGii92l IGBOOI. 
KMC/SSA methods are exact i.e. they are designed to account for every possible discrete reaction event and be- 
come notoriously inefficient for the scenarios when stochastic effects are to be ignored and model can be abstracted as 
system of conventional ordinary differential equations (ODEs) and ambitions goals which require the intensive use of 
Kinetic Monte Carlo methods will also require immense computational resources. That became a main motivation for 
the development of the hybrid simulation methods which will allow adequate levels of description for different parts 
of complex system in hand. 

Recently several attempts were made to bridge the gap between different levels of discretion and in particular 
reduce excessive information content delivered by KMC/SSA methods on a temporal scale. To overcome excessive 
level of resolution on the temporal scale i.e. to perform the temporal coarse-graining r-leaping methods were devel- 
oped IIGilO 1 1 IGP031 IBTB04I IRPYG03I . Leaping methods allow to improve the efficiency of the standard KMC/SSA 
while maintaining acceptable losses in accuracy by approximating the time-inhomogeneous Poisson process counting 
the reaction events via the Gaussian process/diffusion approximation. The key idea here is to take a larger time steps 
and allow for more reactions to take place in that step, but under the proviso that propensity functions do not change 
too much in that interval. When it comes to practical usage of the r-leaping method a robust strategy is necessary 
for deciding when during the coarse of the simulation when to step exactly using KMC/SSA and when to leap ap- 
proximately with diffusion methods |Kur71 GilOl EK86|. When leaping, consideration of the parameter values and 
error-reduction scheme is also necessary. 
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This problem becomes of practical importance because it is often desirable to allow some chemical species to 
be treated as continuous random variables and some to be treated discretely. This is particularly true for the case of 
the transcriptional regulation by transcription factors (TFs). In this case there may be as few as one DNA/TF binding 
site and messenger RNA (mRNA) abundance may be as small as tens of molecules while there may be hundreds 
or even thousand molecules of regulatory proteins per cell. In addition to that, biochemical networks, as probably 
every complex system, show remarkably interesting dynamics due to the presence of different time scales describing 
the evolution of discrete and continuous components. The technical difficulty with implementing the hybrid schemes 
is that standard KMC/Gillespie approach requires constant transition rates between reaction events of the discrete 
species. This may not be the case if some of chemical species are evolving continuously in time or with different time 
scales. 

Several attempts have been made to illustrate the relevance and feasibility of hybrid algorithm, especially for the 
case of highly developed separation of time-scales. Rao and Arkin |RA03 1 consider the case of two groups of species 
evolving with two very different time scales and present stochastic averaging procedure, constructing effective rates of 
the "slow" reaction in the presence of the "fast" species, termed stochastic quasi steady state approximation (QSSA). 

In their work |KMS04| Kiehl et. al. are developing an approach to integrate ODE/SDE schemes with Kinetic 
Monte Carlo schemes. |KMS04| used synchronization algorithm to produce the numerical scheme which bridges 
together asynchronous in nature KMC/SSA algorithms and synchronous ODE integration schemes. 

One should probably agree that following the paths presented in the literature up to date will give solutions which 
are very close to the exact dynamics in the limit of a very small time step, largely different time scales, etc. but verity 
of different implementation approaches presented in the literature to this point indicates the fact that unifying solid 
mathematical description has to be thought to justify the development of the hybrid simulation schemes applicable 
for modeling stochastic phenomena in chemical reaction networks. Not many stochastic algorithms presented in the 
literature so far have gone through the rigorous error control analysis or even search for the fundamental justification 
of the presented methods. 

In this publication we aim to provide not only the framework for the for the hybrid simulation method but also we 

outline it's rigorous mathematical justification including error analysis. 
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Outline of this paper is following. In Section|2]we briefly outline some basic notation and assumption used in 
this article. In the next Section|4]our goal will be to give precise mathematical formulation of the stochastic process 
corresponding to the hybrid representation of the reaction network. We model the stochastic dynamics of the reaction 
network as a dynamics of the composite pair (X, er), governed by the system of the jump-diffusion SDEs. This 
section will be followed immediately by simple example and view on error analysis. Article is ended with discussion 
and outline for the future work. 

2. Basic Notations and Assumptions . Consider the reaction network which represent the dynamics of N well- 
stirred molecular species {Si, i = 1 . . . N} which inter-react through R(r = 1, . . . , R) channels of chemical reactions. 
The current state of the system is completely specified by the A-dimensional vector S = (Si, S2, ■ ■ ■ , Sn) consisting 
of non-negative integer numbers, where each number Si defined to be the current total number of molecules of type 
Si in the system. Formally reaction network can be expressed as a set of the transition rules: 

r = 1 : i/^Si + u^S 2 + . . . + v^i^n ffiSi + ^S 2 + . . . + v N1 S N , (2.1a) 
r = 2 : v£ 2 Si + ^S 2 + . . . + 1^2 Sjv + ^2~ 2 S 2 + . . . + v^ 2 Sn, (2.1b) 

r = R : "irSi + ^ 2 + R S 2 + . . . + v^rSn + U 2R^ + ■ • • + v nr$n (2.1c) 

Transition rates a r (S) : Z N —> R + are the probabilities that reaction event corresponding to the reaction 

r : S — >• S + v r 

will take place in infinitesimal time interval [t, t + dt) somewhere inside the reaction volume. Vectors v r = v~ — 1/+ 
are the stoichiometric Al changes of number of molecules of types Si due to the reaction r. As one can see, each 
reaction channel is specified by two quantities: transition rates a r and stoichiometric changes in molecular composition 
of the system due to each individual event in reaction channel v r . 

In Eqn. ( 12. It k r are transition probabilities for the reaction channel r for any particular configurations of the 
interacting molecules on the right hand side of ( 12. H . The functions h r (S) : 2^ — > Z are defined to be the number 
of distinct reaction combinations of the reacting molecules available in the state S. The overall transition rate can be 
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expressed as: 

a r (S) = k r h r (S), (2.2) 

For the mass-action kinetics transition rates h r are found to be a polynomial functions of variables S. For the simple 
reaction system: 

Si ^ S 2 , 2S 2 h 

statistical weights hi are: ft-i(S) = Si, /i 2 (S) = ^S 2 (S 2 — 1). We refer reader to the sources in lvK81l and (Gil92| 
for more details. 

3. Partition of species and reactions. Our interest in the next step will be to consider different reactions in the 
system in a different manner. We partition the set of molecular species into two groups: "large population group" 
(with appropriate quantities labeled with "X" ) and discrete group (labeled with "a"). Species from the first group 
will be represented by the set of i = 1, . . . , Nx variables taking values on K + : 

Xi , Xi , . . . , Xn x , X.- t G M + 

And present in large copy numbers: Xi ^> 1, Vi The reminder of the N a = N — Nx species represent the discrete 
components of the network: 

<7i,er 2 , • • ■ , &i G ^+ 

Partition of the species into the two separate groups correspond to the partition of the stoichiometric vectors v r into 
two sets : 

which correspond to the changes of species in sets X and a due to the reaction event in the channel r. 
Based on this partition of species we can introduce three groups of reactions: 

(i) First group of reactions which we shall call 72-1, consist of reactions with large combinatorial weights, 
i.e. which satisfy the condition h r (-) ^> 1. We state that this condition is sufficient for the use of the diffusion 
approximation (see Appendix SeclAl 
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(ii) Second group , IZ2 consists of reactions which do not satisfy the constraint h r ^> 1 but change the species 
of the group X, i.e. those reactions for which ^ 

(iii) Finally the third group, which we call consists of reactions for which vf, ^ 
Note that, those subset of the reactions are not necessary independent,!. e. 

Tlx n n 3 ^ 0, n 2 n n 3 + 

Without the loss of generality we shall assume that reactions from the set 7?-2 U 7^3 (set of "jump" reactions) are formed 
by reactions in Eqn. \2.\ \ with indices 1, . . . , Rd, Rd < R- The rest, R — Rd are assumed to satisfy the condition of 
having large combinatorial weight. 

4. Formulation of stochastic dynamics: jump-diffusion models. Chemical Master equation lvK81l lGar85l 
is usually employed to describe the stochastic dynamics of reaction networks of the type ( 12. 11 1. It is notoriously hard 
to solve analytically even for the simple systems | SR95 1 and, as we mentioned before, one usually has to resort to so 
called Kinetic Monte Carlo methods IGU77I IBML75I to obtain numerical solution. In this article we shall abandon 
the description based on the Chemical Master Equation(CME) and resort to the explicit treatment of the stochastic 
trajectories of the process {St} as a solution of the system of stochastic differential equations (SDEs) providing the 
recepie for the construction of numerical schemes. 

Poisson counter process is a simple but important process and below we extend its utility by combining with ideas 
from differential equations. We briefly describe this approach to pave our way to specific applications. 

To describe the structure of the stochastic jump processes let us introduce in the rather traditional way a stochastic 
basis (f2, J 7 , P) supplied with the filtration J-~t>o large enough to include i?-dimensional unit Poisson processes. We 
recall that state variables S t are defined on this probability space. On this space transition rates Eqn. ( 12. 21 define the 
set of time non-homogeneous counting processes N r (t)\ GS79 1 : 

E(N r (t + At) -N r (t)\T t ) = a r (S t _)At + 0(At 2 ), (4.1a) 

Notation <_ is used to underline the existence of jumps in the process St', t- = lim^o^ — S). The above equation 

states that differentials of processes N r (-), dN r = N r (t+dt)—N r (t) depend upon local intensities a r (St_ ). Notation 

t- is used to underline the existence of jumps in the process S t : t_ = lima^o(^ — §)■ Counting process N r (-) can be 
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also represented as a time changed unit-rate Poisson processes II r (-), r = 1, . . . , R [EK86|: 

N r (t) = n r (^J a r {S s _)ds^j (4.1b) 

In turn, counting processes N r (-) represent the stochastic dynamics of the Markov process St through the "mass- 
balance" equation: 

R 

S t = S + Y,»rN r (t), (4.2) 

r=l 

Solution of Eqn. ( 14.21 is a constant on the interval where all N r (•) are constants and jumps by v r whenever the reaction 
event takes place in the reaction channel r at the moment of time t: 

S t = S t _+i/ r (4.3) 

Using the partitioning which was introduced earlier, stochastic dynamics of the state vector (X, er) can be presented 
as following system of SDEs: 

R 



X t = X +5>*iV r (i), (4.4a) 

r=l 
R 

<r t = <ro + X) l W(t), (4.4b) 



where integer vectors and v° represent the change of the components X and a due to the reaction event in the 
channel r. In general, propensity functions a r ( ) depend on both X and er. 

Large combinatorial weights h r {-) 3> 1 of some reactions allow us to follow the diffusion approximation |Gil01 
EK86|, i.e. by using the weak convergence methods we express differentials of the time-inhomogeneous Poisson 
processes in Eqn. ( 14. 2> corresponding the group 1Z\ as a Gaussian processes with the drifts a r (-) = a r (X t , er t ): 



dN r (-) -» a r (-)dt + y/a r (-)dW r (t), (4.5) 

where W r (t),r £ IZi are Brownian motions independent of each other and X.t=o ■ This approximation of the stochastic 
dynamics allows us to make a simplification in which we can consider the components (X, er) as defined on a hybrid 
space M^ x x ll^" ; X takes values on the continuum state space, , while er still takes values on discrete state 



spaceZ + ". This brings us to the following formulation of the stochastic dynamics of the pair (X, <x) 



dX t = i^Or(X t _,o- t _)dt+ <ay 2 (Xt_,*r t JdW r (t) 



dif fusion 



- ^ d/V r (a r (X t _,o- t J), (4.6a) 

r£K 2 



jump 



da t = J2 <diV r (a r (X t _,o- t J), (4.6b) 



where we explicitly demonstrated the dependence of the differential of the counting processes dN r (-) on local inten- 
sities. Note that according to the presence of the jump component term in Equation (|^^aj X t ^ X t and X t is a 
jump-diffusion process in . 

Processes with state-dependent intensities of jump are notoriously hard to model. One complications which makes 
the study of the SDEs of the type Eqn. J4.6i hard is the nature of the noise source driving the dynamics. In contrast to 
very well studied Wiener or Poisson driven SDEs | KP92 1 is that intensity of the noise source is stochastic in nature. 
We shall take a different representation of the noise source. Namely we use the technique of the Poisson random 
measures. Our discussion of the random measures will be brief and rather informal; for mathematical background see 
monograph by Jacod and Shiryaev | JS87 1 or Cox and Isham |CI80|. 

Poisson random measure fi(ui,dt x dz), ui G ft defines a sequence of points and marks {(rj, Zi)} with the simple 
interpretation that the mark Zj arrives at time r^; Tj take values on K + . Marks Zj take values in a general space 
E , which in the context of the problem considered here can be taken as a subset of the Euclidian space, E C KL 
Deterministic intensity of the Poisson process n(dt x dz) is m(dz) dt with m{dz) being a measure on E. Arrival 
times follow a Poisson process with deterministic intensity dt and marks a i.i.d distributed in the interval dz. 

Important procedure which can be used to simulate the intensity dependent process Eqn. \2.\\ is based on in- 
troduction of the partition on the marked space. At every point of the state space S = (X, er) we construct the 
set of disjoint intervals based on transition rate of "jump" reactions, i.e. reactions from the subset 7?-2 U 7?.3, i.e. 
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r = l,...,R d : 

A r (S) = [A r _ 1 (S),A r (S)), (4.7a) 

r 

A r = 0, A r (S) = a r(S) (4.7b) 

r' = l 

where index r runs over reactions which either do not satisfy condition h r ^> 1 of for which vector v a r . is non-zero i.e. 
they change discrete components. Thus in general length of the interval A r (S) is a r (S). Now define a set of functions 

c r : xZf xl^ {0, 1}, r = l...R d : 

!1, if z e A r (X,er), 
(4.8) 
0, otherwise 

Suppose now that the space (f2, J 7 , J~t, P) supports a Poisson random measure /J,(u>, ds x dz) = /i(dtxdz), w £ f2 in 
R+ x E, £cR with a compensator dt dz and a set of Wiener processes W r (-), r S 72. i. Then time non-homogenous 
counting process N r (t) — J Q dN r (a r (X s _ , cr s )) can be represented as BGS79IIPro83l : 

N r = J J Cr(X s _,a s _,z)fj,(ds x dz) (4.9) 

and (R+ c x Z^ d ) valued process (X t , cr) can be viewed as a solution of coupled Ito-Skorohod SDEs driven by the 
Poisson random measures |Pro83 GS79|: 

dX t = J2 ^a r (X t _,o- t _)di+ <a 1 r /2 (X t _,<T t _)dW r (t)+ 

reKi reTZ ± 

+ X] / K c r(X- t _,(T t _,z)n(dt x dz), (4.10a) 
dcr t = ^ / L>°c r {X t _,(T t _,z)fj,(dt x dz) (4.10b) 

Function c r in Eqn. J4. lQi transforms arrived marks into magnitude of the jump of S due to appropriate reaction 
channel. 

For the practical insight on the dynamics described above, we will introduce reference Poisson process. 

Since the length of each interval A r (S) is a r (S) and bounded function for every r = 1 . . . R then it is follows 
that the length of the interval A(S) (we denote it |A(S)|) is bounded function of S = (X, cr). Therefore, it has a 
maximum at some point S* = (X*, cr*): 

|A(S)| < |A(S*)| (4.11) 
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Let A max = A(S* ) | denote the maximum length of the joint intervals in the marked space. Now consider a standard 
Poisson process with intensity A max . Let r^,n — 1,2,... denote the jump times of this process (sampled, perhaps, as 
the partial sums of the exponentially distributed random variables with mean A"^). Let A* = A(S*) be the marked 
space and {z n } be the sequence of i.i.d. distributed random variables with uniform distribution on A*, independent 
on N t . In this case we can represent a Poisson random measure with intensity dz dt is related to the marked point 
process (t^, z„)„>o, constituting the sample path of the process. Namely for each A £ A*: 

JV A ((0,t], A) = lrA< t l* n eA, (4.12a) 

n>\ 

E(AT A ((0, t],A)) = A max tF(z n € A) = t(i{A) (4.12b) 

This representation (sometimes referred to as thinning of the Poisson measure |CI80|) is very useful in practical 
numerical applications. In the next section we will use the reference process N\ to construct the discretization scheme 
taking A — A r (-) ( 14. 7> for various r £ 72.2,3' 

5. Implementation of Numerical Scheme. In this section we outline basic principles for the construction of 
the numerical scheme for the solution of the system given by Eqn. j4.10l i. To begin with, we should introduce the 
appropriate discretization of the interval [0, T] on which dynamics of the system has to be investigated. Let us denote 
by [0, T]h t M = {0 = to < t\ < t% < . . . < tM+i = T}, perhaps the usual equidistant, discretization of the time 
interval , i.e. h = T/M. 

Suppose now that {0, r^, , • ■ ■ }, T 4 A £ [0, T] are jump times of the reference Poisson process l4.12al : then we 
consider new time discretization which is a merger of the points from [0, T]^^ and jumps of the reference Poisson 
process Pi. 12al 

It is important to note that jump-times of the Poisson measure Eqn. ( 14. 12ai can be modeled without discretization 
error and with little prior information; one has only make sure that transition rates from the groups 7^2,3 are bounded, 
i.e. are less then chosen intensity of the reference Poisson measure. Armed with the time discretization scheme 
assembled from equidistant discretization and jump-times of the reference process we introduce the following simple 
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jump-diffusion numerical scheme: 



Xj,ti — Xj,ti-! + ^ ^a r (Xt i _ 1 ,o- 4i _ 1 )(t i — tj_i)+ 

r£Ri 

+ i^oJ/^X^.^JCWr^O-W;^-!)), (5.1a) 

X j)tj = Xr t . + vfr j £rptu,*U-n z )Kdt x dz), (5.1b) 

*i,t< = Vj,U-i + ^2 V jr / ^feti^U-iiZ) M d * x dz ) ( 5 - lc ) 



which recursively determines the values of the discretized processes X = (Xi, . . . , Xjv ) and <r = . . . , <7jv d ) 
at points ti, starting from values Xq, <tq. Increments of the Winer processes W r (-) are zero-mean Gaussian random 
variables: W r (ti) — W r (tj_i) oc Af(Q,ti — £i_i). According to the i5.l\ . between the jumps dynamics of the 
component X is purely diffusive. Integral in Eqn. (15 . 1 bi is evaluated at the single point Zi which belongs to the 
uniform distribution U([0, A max \), i.e. jump in channel r is accepted iff generated random variable z ti £ ^([0, 1]), 
generated on the step ti satisfies the condition: 

ZtAmax £ KiiXJ^&U-J), (5.2) 



Depending on the relative ratio between characteristic time-scales of diffusion and discrete species different nu- 
merical schemes must be considered to achieve substantial speed-up of the simulation. Here we consider only the 
case of fast diffusion modes and slow discrete modes. In this case every interval [r A , r^ +1 ) should be partitioned onto 
smaller intervals of the discrete grid on which numerical scheme (15. 1> for diffusion approximation is used. Opposite 
case of the fast switching in discrete component er will be considered elsewhere. 



6. Example. We now wish do demonstrate application of the above numerical scheme. In this section we con- 
sider the simple example, which nevertheless, is interesting enough for the purposes of the demonstration of the 
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technique discussed in this paper. Consider the following model of the reaction network: 

Si^Sz, (6.1a) 

S 2 h S l5 (6.1b) 

Si^Si+nS 3 , (6.1c) 

S 2 + S 3 h S 2 + S 4) (6. Id) 

S 4 ^ (6.1e) 

with continuous and discrete components <x = (Si, S2) = (01,02) an d X = (53, S4) = (Xi,X 2 ) respectively. The 
functions a r and vectors v r for this system are: 

ai(X,<r)=fti<7i, v\ = (-1,1,0,0), (6.2a) 

o 2 (X, 0) = k 2 <7 2 , v\ = (1, -1, 0, 0), (6.2b) 

a 3 (X,cr) = fc 3 cri, i/f = (0,0,n,0), (6.2c) 

a 4 (X, cr) = k 4 02Xi, v\ = (0, 0, -1, 1), (6.2d) 

a 5 (X, <t) = fc 5 ^2, i/f = (0, 0, 0, -1) (6.2e) 

We have used the following set of kinetic parameters^ = 0.50 , k 2 — 0.50 , fc 3 = 1.00 , fc 4 = 0.10 , fc 5 = 0.01 , 
Xi(0) = 1000, X 2 (0) = 200 and n = 5 1 . 

Analysis of the propensity functions (/i 4;5 3> 1) shows the following partition of the original reaction set 1Z = 
{1.....5}: 

^i = {4,5}, (6.3) 

^ 2 = {1,2}, (6.4) 

Tl 3 = {3} (6.5) 

i.e. reactions {4, 5} correspond to the diffusion modes. We also constraint ourselves in this example considering the 



'One can consider this as a simple model of transcriptional regulation, where presence/absence of the transcription factor (Si G {0, 1}) leads 

to the bursts of the transcription of protein S3 with n proteins per transcription event. 
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case: 

ai + a 2 + a 3 = A max < a 4:5 (6.6) 

Propagation of the diffusion mode was performed via Euler-Maruyama scheme |BTB04|,[KP92 1 with set of different 
time-steps: h = 0.1, . . . , 2.0. We performed 10 4 Monte Carlo runs with both standard KMC/SSA and hybrid scheme 
outlined above to obtain the distributions of the number of molecular species (S3, S4) — (Xx^X^) at time-slice 
T = 2 x 10 3 . 

Figures lC"TlC"2l show that distributions of the species S 4 , S 5 , P(S 3 ,T) = P(X 1 ,T), P(S 4 ,T) = P(X 2 ,T) 
obtained with the hybrid scheme is almost indistinguishable from the results of the obtained by following complete 
KMC/SSA-based approach. We also report an observed speed-up in terms of the ratio of the CPU times I, KMC , 

J- hybrid 

Fig. (IC.3> as a function of time step h used in numerical integration of diffusion modes. 

7. Discussion. Advances in field of molecular biology outlined new directions of research in stochastic chem- 
ical kinetics. We now need to develop the software tools and mathematical approaches to integrate models from 
micro-scales to macro-scales in a seamless fashion. Such multi-scale models are essential if we are to produce quanti- 
tative, predictive models of complex biological behaviors. An important step on this path is implementation of hybrid 
simulation methods, capable to account for heterogeneity of properties of interacting components. 

We have outlined rigorous framework for development of hybrid simulation schemes based on the path sample 
representation of reaction dynamics rather then on CME based approach. For the best of our knowledge this type of 
analysis has not being performed before. 

It was not the primal goal of this publication to justify the partitioning of the species/reactions into the different 
groups. One should probably have some prior knowledge when partitioning the species into different types. Choice of 
the reactions ("diffusive" or "jump") is based on functional law of large numbers and error generated by this step can 
probably be controlled rigorously with some insights on this subject outlined in SectionlAl 

Numerical solution for the propagation of the diffusion modes were outlined via simple constant step explicit Euler 
scheme. Typical explicit or implicit scheme |RPYG03| based on stochastic Taylor approximation |KP92 |,[Mil95 1 
might run into problems , since they do not conserve the non-negativity of the numerical solution, i.e. i.e. it cannot 
guarantee Xj^ G K+ almost surely. We point out here that use of balanced-implicit stochastic schemes (BIM), 
|MPS98| look promising when it comes to the construction of the numerical solution with the property of "almost 
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sure" positivity. BIM scheme known to have the same order of convergence as the Euler-Maruyama stochastic scheme 
, namely , the error OQi 1 ' 2 ) for approximations of the individual trajectories and O(h) for moments. 

Adaptive time-step control systems also look promising as a direction for the research aimed for achieving the 
simulation speed up|BTB04|. 

8. ACKNOWLEDGMENTS. Author would like to acknowledge valuable discussions and useful comments 
from Prof. C. Myers, Prof. D. Steinsaltz, K. Ericksson and Prof. A. P. Arkin. 

Appendix A. Diffusion limit of the Poisson process with state-dependent intensity . 

Introducing the set of independent Poisson processes with rate k r differential of the counting process dN r at some 
state St can be expressed as: 

n=h,.(S) 

dN r = dU k} = k r h r (S)dt + dM r , (A.l) 

i=l 

n=h r (S) 

dM r = dnjjj = <mf r - k r dt, E(dnjJ|.Ft) = o (A.2) 

i=i 

Using the functional analog of the law of the large numbers one can conclude: 

dM r = y/k r h r dW r + de(t), (A3) 

where E(||rfe r (t)|| \\F t ) oc h^ 1 (S) < 1 as h r > 1 (see IIEK86I IAK95 1). Note that this condition is different from 
previously reported in the literature | GP03 1 . 

Appendix B. Convergence analysis of the Hybrid Approximation. 

Below we demonstrate convergence properties of the numerical scheme J5.ll5.lbf if propensity functions a r (-) 
satisfy regular Lipschitz continuity condition: 

MX,er) -a r (X',cr')| 2 + | a y 2 (X,cr) - <4 /2 (X', a')\ 2 < 

< K(\o- - o-'\ 2 + \X-X.'\ 2 ), Vr eTZ, (B.l) 

where |X|,|er| are the Euclidian norms of the vectors X,er respectively and K is some constant. 

Since in scheme l4.6bl generation of the jump times of the counting process N\(t) is not linked to the dynamics 
of the state space pair (X, cr) we can consider first the error of the approximation of the scheme d4.6bi for a given 
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sequence of the jump times of the Poisson process (14. 12ab N\ : 

< < r 2 A < . . . < < T, N A (T) = N, 

with random number N following the Poisson distribution: 

(A T) N 

P(N A (T) = N) = [ > exp(-A maa T) (B.2) 
On the interval [t^, t£ +1 ) we represent the numerical solution corresponding to the diffusion mode as: 

X t = X t a + / v r a r (X. s , <T T ^)ds + v r al^ 2 (s, X s , & T A.)dW r , (B.3) 

X f = X ffc , t G [t k ,t k+ i) C [Tk>T~k+i), X t=0 = X , <T t=0 = er . (B.4) 

On each interval between two jump times of the reference process [t a , t^ +1 ) we perform our analysis in a rather 
traditional way, starting from estimating the expectation of the norms |X t — X t | 2 , \a t — <x t | 2 conditional on filtration 
generated up to jump time r A : T t k : 

e x ([r„,T n+1 )) = sup E(|X t -X t | 2 (B.5a) 

e CT ([Tn)7"n+l)) = sup E(|er t - <x f | 2 \T t a) (B.5b) 

*£[^,r A +1 ) 

Then, evidently, total error on the interval [0, T] for each of the components, X and er can be found as: 

e x ([0,T)) = maxe x ([r n ,r„ +1 )), (B.6) 

n 

e CT ([0,T)) = maxe CT ([r„,T n+1 )) (B.7) 

n 

Let us turn to the estimation of JB.5a> on each interval [t a , t„ +1 ) based on the formulation ( IB.3I >. For the component 
X for every t G [r„ , t^ +1 ) one has: 



ex([0,i] 



sup E([X Tn -X Tri + / V" f r (a r (X„cr„) - a r (X u ,o-„))^ 

S&[0,t] JO „^T,_ 



5] ^(ay 2 (X uCTu ) - a y 2 (X u ,«T„))dW r ( U )] 2 |^ T A) < (B.8) 



rERi 



<|X r „-X T J 2 + sup E(|Ai(s)| 2 |Jva) + sup E(\A 2 (s)\ 2 \F t a) (B.9) 

s£[0,t] s£[0,t] 
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where: 

Ai(s)= / i/ r (a r (X M ,(7 T J - a r (X u ,& Tn ))du, 

and 

A 2 (s) = [ Yl Mal /2 (Xu,<T T J - a y 2 (± u ,& T J)dW r (u) 
By Couchy's inequality bound for the second term in lB.9l can be presented as: 

sup E(|Ai(s)| 2 |jr rA ) < (t-r£) sup E( [ ( V i/ r (o r (X u , oyj - a r (X u , <x r „))) 2 <H^rA), (B.10) 

s6[0,t] sG[0,t] Jr,A 



At the same time, by Doob's inequality for martingales, third term can be estimated as: 

3 E(\A 2 {s)\ 2 \T t a) <4E(| J 4 2 (t)| 2 |^ rA ) = 4E / 

se[o 



sup E(\A 2 {s)\ 2 \T t a) <4E(| J 4 2 (t)| 2 |^ rA ) = 4E f V (a*/ 2 (X u , o- T „) - ay 2 (X n , a T J) 2 du (B.ll) 



which together with Lipschitz condition {Bli and Gronwall lemma gives the bound for the error of the component X 
on interval [r A , t], ex([r A , t]): 

ex([T^,t])<K f ex([T^,s})ds + KTE{\± t a -X t a\ 2 + KTE{\ct t a - & t a\ 2 ), (B.l2a) 

J T A n n n 

£x([Tn,t}) <KT[E(\X t a - X t a\ 2 + KTE(\ct t a - <t t a| 2 )] exp(K(t - r A )) (B.l2b) 
One can see that additional error due to the possible jump at time r^ +1 terms is also bounded: 

/ l^( X rA +1 _^rA +1 _,^) -c,.(X r A +i _,(T r A +i _,z)| 2 dz < K(\X t a +i _ - X t a +i J 2 + |<x t a +i _ - <x r A +1 J 2 ) , Vr G TZ 2 . 

J IS 

(B.13) 

Assuming that one is using scheme with p-order of strong convergence (p = 1/2 for Euler scheme) | Mil95l : 

sup E{\X t -X t \ 2 + \a t -& t \ 2 \T Q )<K(l + \a a \ 2 + \X a \ 2 )h 2 ^ (B.14) 

te[o,TA) 

and using this inequality together with ( IB.12b> . JB. 1 3i one can come to the conclusion that the overall strong error of 
the scheme is 

ex([0,T))<C({r I f},r)/i 2 P, (B.15) 
e a ([0,T))<C({4},T)h 2p (B.16) 
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where constant C depends on the sequence of jump times {r^} and length of the interval T. Hence after averaging 
over all possible jump-times overall error has a strong asymptotic of the numerical scheme used to propagate the 
diffusion modes, i.e. e x ([0, T)), e ff ([0, T)) oc h 2 P. 
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